Nonlinear rupture of thin liquid films on solid surfaces 
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In this Letter we investigate the rupture instability of thin hquid films by means of a bifurcation 
analysis in the vicinity of the short-scale instability threshold. The rupture time estimate obtained 
in closed form as a function of the relevant dimensionless groups is in striking agreement with 
the results of the numerical simulations of the original nonlinear evolution equations. This suggests 
that the weakly nonlinear theory captures the adequate physics of the instability. When antagonistic 
(attractive/repulsive) molecular forces are considered, nonlinear saturation of the instability becomes 
possible. We show that the stability boundaries are determined by the van der Waals potential alone. 
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PACS numbers: 47.20.-k, 68.15+e 



It is well known that a liquid film on a planar solid sur- 
face may become unstable due to long-range molecular 
forces. The forces originating from van der Waals attrac- 
tions m accelerate thinning in regions of film depression 
leading to film rupture and "spinodal dewetting" |2| . On 
the other hand, electrical double layers on the solid sur- 
face may give rise to intermolecular repulsions stabilizing 
thin films against rupture |3|. 

In recent years, much effort has been put into theo- 
retical modelling the dewetting phenomena |j, S B S 
H, M, O H 01 . A nonlinear theory of the film evolu- 
tion based on the long-wave nature of the response was 
first posed in Ref. y. This approach, which has already 
been considered for different situations |5|, yields nonlin- 
ear partial differential equations that describes the evolu- 
tion of the interface shape, surfactant concentration, etc. 
Linear stability analysis is routinely applied to predict 
the onset of the instability and the characteristic wave- 
length, but the rupture time estimate obtained from the 
linear theory turns out to be rather poor: it underesti- 
mates the rupture time due to highly nonlinear nature 
of response. The most common and straightforward ap- 
proach is to solve the evolution equations numerically [g , 
113 > El El HI 16]- The obvious disadvantage of the 
numerical simulation is that for a complex problem that 
involves many parameters, full parametric study of the 
rupture is quite elaborate. 

A bifurcation technique was first applied in [g to ar- 
rive at the nonlinear estimate for the rupture time in the 
vicinity of a steady bifurcation point. It was demon- 
strated that nonlinear terms owing to van der Waals 
attractions contribute to rapid acceleration of the rup- 
ture beyond the linear regime. Analysis of the nonlinear 
evolution of small disturbances leads to a dynamic Lan- 
dau equation for the perturbation amplitude. The closed 
form solution of the amplitude equation provides a time 
for "blowup" of the initial small-amplitude disturbance 
that was proposed to be a good estimate of the nonlinear 



rupture time. The approach has never been given enough 
attention perhaps because the analysis involves rather te- 
dious algebra and can only be done "by hand" for some 
simple cases. It has been demonstrated in |9| that the 
derivation of the amplitude equation can be automatized 
by using a previously developed symbolic algorithm for 
bifurcation analysis [l(]j. Although the closed form non- 
linear estimate of the rupture time of the thin film in 
presence of insoluble surfactant was derived in |9|, the 
lack of parametric study of the problem by simulations 
didn't allow a proper comparison of the two approaches. 
Recently, an extensive numerical study of the thin film 
rupture driven by van der Waals attractive forces in the 
presence of insoluble surfactant and hydrodynamic slip 
was reported in |ll| . We have developed a generalized 
theory of thin film rupture for an arbitrary intermolecu- 
lar potential; further, we compare the rupture time esti- 
mate from our theory with the results of simulations by 
|ll| for the purely attractive potential and come up with 
some predictions for the competing (attractive/repulsive) 
potential. 

We consider a model describing the evolution of a thin 
liquid film a solid substrate subject to a van der Waals 
force in the presence of a slip and insoluble surfactant. 
The dimensionless film thickness h and surfactant con- 
centration F are governed by a system of coupled evo- 
lution equations derived in the long-wave approximation 

Ell, 



ht 
with 



MT^h [-+f3] - J'Ji' 



-/3 



i XX 

~P' 



MTT,{p + h)-TT,h[-+p 



(1) 
,(2) 



T{x,t) = -^{h) +Chxx 
and where A1 is a Marangoni number, P is a Peclet num- 



ber, /3 is a Navier slip coefficient, C is a surface tension 
parameter and </? = (dAG/dh) is the van der Waals po- 
tential (all dimensionless) . For nonslippingfilms {(3 = 0) 
equations (|1I2|I are equivalent to those in J3- 

The linear stability analysis of the uniform station- 
ary state Uq = {hojTo} results in a critical value of 
the wavenumber corresponding to a stationary bifurca- 
tion point kc = {—ip'{hQ)/Cy/^. Following a standard 
procedure we choose C as a bifurcation parameter. In 
a bounded domain, < x < C, the basic solution 
Uq changes stability (becomes spinodally unstable) with 
k = kc and when C < Cc — — (/?'(/io)>C^/47r^, where Cc 
correspond to a steady bifurcation point. 

To investigate the nonlinear problem in the vicinity of 
the bifurcation point we expand the bifurcation param- 
eter as C = Cc + e^C2 + ■■■, where e is a small critical- 
ity, introduce a slow time scale suggested by the linear 
theory, t = e^i, and seek the solution in power series 
of e as u = Uq -I- eui -I- .... Substitution of this expan- 
sion into the system (|1I2|) to the first order in e yields 
ui = (yl(T)e''^<=^ -I- cc.) U, where U = {1,0} is a solu- 
tion of the linearized zero-eigenvalue problem. The com- 
plex amplitude A satisfies the dynamic Landau equation 
which is determined to 0{e^) of the perturbation theory: 
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where the linear coefficient a and the Landau coefficient 
K are given by 



a = -S kfC2, 



respectively, and 
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The closed form solution of the amplitude equation 
can be easily obtained given the initial value of the am- 
plitude Aq = A{Q). The "blowup" time, corresponding 
to the infinite growth of the amplitude A and providing 
a nonlinear estimate of the rupture time (in the original 
time scale), in the vicinity of the bifurcation point, as 
C2 ^ is determined solely by the cubic coefficient 
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where Aq is now 0{e). An important observation is that 
despite the complicated nature of the original evolution 
equations, the structure of the cubic coefficient k in Q 
is very simple: terms in brackets contain only deriva- 
tives of the intermolecular potential (p and a factor 5 /Cc 
incorporates the dependence on the rest of parameters. 
Since 5 > 0, it is readily seen from Q that for purely at- 
tractive potential k is always positive and the rupture is 
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FIG. 1; Variation of rupture time vs. /3 with Ao = 0.106, 
ho = 1, C = 1, M = I, V = 100, To = 0.5, c = 3 and ^ = 1 

( ), A — 2 ( ). The inset shows analogous results of 

the numerical simulations Il3. 



inevitable. For the most commonly encountered attrac- 
tive potential 4> = A/h'^, with c = 3,4 (unretarded and 
retarded case, respectively) we calculate the rupture time 
from l^l and compare to the results of numerical simula- 
tions of the original evolution equations (|1I2|I reported in 
[ll|. A = A^/Qnpi'^h'^''^ is the scaled Hamaker constant 
with h^ being the mean film thickness (here and there- 
after dimensional quantities are marked with ^). Typical 
evolution of the film thickness in numerical simulations 
shows an accelerated thinning of the film in the depressed 
region due to the initial disturbance at some unstable 
wavelength A > 27r/fcc with and subsequent film rupture. 
This behavior suggests that the acceleration of the film 
rupture is due to a nonlinear self-coupling of the pertur- 
bation beyond the linear regime and therefore the cubic 
nonlinearity in the amplitude equation ||2Jl should provide 
an accurate description of the nonlinear rupture. Since 
^ is formally valid in the close vicinity of the instabil- 
ity threshold, it is considered as a nonlinear approxima- 
tion for the rupture time, tmp, far from the instability 
threshold (as in [llj) with one adjustable parameter, Aq. 

Figures HEl show the comparison between the theory 
with Ao = 0.106 and the simulations [llj. It is read- 
ily seen that there is an excellent quantitative agreement 
with the numerical estimate of tmp as a function of dif- 
ferent dimensionless parameters. An obvious advantage 
of the present analysis is that the closed form expression 
for the nonlinear rupture time as a function of different 
parameters of the problem is available in a closed form 
expression for a general van der Waals potential. 

When the antagonistic attractive/repusive intermolec- 
ular interactions are present, nonlinear saturation of the 
rupture instability is possible as k may change sign (su- 
percritical bifurcation). Let us consider the general rep- 
resentation of the antagonistic van der Waals potential 
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FIG. 2: Variation of rupture time with /3 for A = 1 and 
the same values of other parameters as in Fig. except we 

vary M: M = I ( ), M = 100 { ). The inset shows 

analogous results of the numerical simulations [l^ . 
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FIG. 3: Variation of rupture time with (3 for A = 1 and 
the same values of other parameters as in Fig. Q except now 

we vary C: C = 1 ( ); C — 0.1 ( ). The inset shows 

analogous results of the numerical simulations [l^] . 



miii 






■ {Sp/lp)e^p{-h/lp), (6) 

where A is defined as before, B = B^, / pf"^ hf~'^ , Sp — 
Sp^hl/pu"^ and Ip = Ip^/h^. It follows from Q that 
the stationary nonruptured solution with amplitude A = 
(— o/k)-^/^ is stable if 



d7i + if 72 < 



EE 

ij; i+j<2 



^tjllli < , 



(7) 



where the first inequality is imposed by the linear theory, 
71 = Bhl-'^/A, 72 = {Sp/lp)hle-"/A, H - ho/lp and 
Uij are some polynomial functions of c, d and H only. 

When the interplay between algebraic potentials is 
considered (72 = 0), the nonlinear stability region is 




FIG. 4: Stability diagram for an antagonistic exponen- 
tial/algebraic potential (71 = 0) in plane of parameters 
72 = {Sp/lp)h'i,e-'^°''"P /A and H = ho/lp. The dashed curve 
corresponds to marginal stability boundary with c = 3. The 
regions of stationary nonruptured states are shown in gray. 
Colored curves show the dependence 72 vs. H for varying 
hydrophobicity of the substrate from hydrophilic (blue curve) 
to hydrophobic (red curve) using parameters from 15]. 
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FIG. 5: Stability diagram for a general van der Waals poten- 
tial Q . The regions right to the dashed curves are spinodally 
unstable and the regions between the solid curves correspond 
to stationary nonruptured state. Black curves correspond to 
A* = 3.0 X 10~^° J, B, = 5.04 X 10"^^ J m** and l-p^ = 0.6 
nm (nonlinear stability regions are shown in gray). The color 
curves correspond to the same values of the parameters ex- 
cept /p^ — 0.4 nm (blue curves) and A = 1.4 x 10~^° J (red 
curves) . 



defined by 71 alone. For exponents (c, d) = (3,4) (re- 
pulsive retarded van der Waals force) the film is stable 
if 0.51 < 71 < 0.73. For the exponents (3, 9) (short- 
range Born repulsion) the stability window is shifted 
to lower values of 71 and the film is stable whenever 
0.066 < 71 < 0.30. For instance, using the values of 
the Hamaker constants measured for a polystyrene film 
on oxidized Si wafers with A* = 2.2 x lO^^*' J and 
Bt: = 5.04 X 10"^'°^ J m^ JI3 the nonlinear analysis pre- 
dicts a stable film thickness of /iq « 2 nm (71 = 0.068) 



while the hnear theory results in 1.56 nm and the equi- 
librium thickness determined from the minimum of AG 
is only 1.3 nm [l3j. For thicker films 71 is rapidly de- 
creasing as h^^, e.g. 71 « 0.00125 for ho = 3.9 nm. 
Although, a qua litative difference in the morphology of 
dewetting in [13 was observed for thicker films (3.9 nm 
vs. 4.9 nm) we speculate that the qualitative change is 
due to approach to the nonlinear stability boundary as 
/i, approaches the thickness of 2 nm. 

When the interplay between algebraic and exponential 
potentials is considered (71 = 0), the nonlinear stabil- 
ity diagram can be defined in terms of 72 and H as in 
Fig. 01 The dashed lines correspond to the boundary 
of linear stability, while the regions of nonlinear stability 
corresponding to k < are shown in gray (iii and iv). 
Note that, if the long-range attraction is combined with 
the shorter-range repulsion {A > 0, Sp > 0), the spin- 
odally unstable region is under the dashed curve (region 
ii in Fig. EJ, and nonlinear saturation occurs for H > 4.5 
and small values of 72 (region iv). In the opposite case 
(^ < 0, 5p < 0) the spinodally unstable region lies above 
the dashed curve (region i) and stabilization occurs for 
thinner films, if < 4, and moderate values of 72 (re- 
gion iii). For instance, for aqueous films on Si substrates 



with c = 3, A^, 



-1.41 X 10 



-20 



J and L 



0.6 nm 



115] we plot 72 vs. H in Fig. ^ (color curves) for dif- 
ferent values of Sp^ varying from —0.61 mJ/m^ (blue) to 
—8.5 mJ/r n^f red) due to increasing hydrophobicity of the 
substrate [l5j. It is evident from Fig. 0]that the emer- 
gence of stationary nonruptured ultrathin films is pos- 
sible on non-hydrophilic substrates as the color curves 
cross region (iii), while on hydrophilic substrates (the 
blue curve) the film of any thickness is stable, in accord 
with 0- 

More interesting behavior is anticipated for competing 
short-range algebraic and exponential potentials. In this 
case 71 , 72 7^ and as they both vary with ho we chose 
to depict the stability diagram in terms of dimensional 
quantities, Sp^ and /i, as in Fig. [S] It is evident that 
stabilization is possible for a wide range of film thick- 
nesses, /i*. When the magnitude of the exponential re- 
pulsion is small, the steady nonruptured state is only 
possible for ultrathin films; for moderate values of Sp^ 
the band of stable solutions widens. For instance, when 
Sp^ — 1.1 mJ/m^, with parameters corresponding to the 
red curve in Fig. nonlinear theory predicts that the 
film is stable below a thickness of ~ 4 nm, whereas linear 
stability provides a value of ~ 2.2 nm. 

Finally, we consider van der Waals interactions of the 
polystyrene films with SiO coating on Si substrate 



m 



dimensional form) 
Asi = -1-3 X 10-19 J 



V'* = — 



-■A,siO 



with 



■A,sio 

20 J and B, = 



hi 



5.04 X 10-^5 J m^ hi 



6ii-(ft,.+d,)3 

A^SiO = 2.2 X 10 

In this case the stability diagram 



is a narrow stability window for ultrathin films up to 
~ 2 nm while its width is insensitive to the variation in 
SiO coating thickness. The stabilization for thicker films 
does not materialize similar to the previously discussed 
case of exponents (3, 9) without coating. 

In conclusion, we have developed a nonlinear theory 
for the rupture of a thin liquid film subject to a general 
van der Waals potential. The comparison between the 
prediction of the weakly nonlinear analysis and the nu- 
merical results is provided for the first time; it is demon- 
strated that there is an excellent quantitative agreement 
between the nonlinear rupture time estimate from our 
theory and the numerical estimate. When an antagonis- 
tic potential is considered, the saturation of the rupture 
instability beyond linear regime is possible, while the sta- 
bility boundary is determined solely by the intermolec- 
ular potential. The results concerning the existence of 
steady nonruptured states should be accessible via nu- 
merical simulations. 



can be depicted in terms of the film thickness, /i*, and 
the SiO coating thickness, d, (not shown). Again, there 
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